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Abstract. The recent years have witnessed a growing interest for covariant 
Lyapunov vectors (CLVs) which span local intrinsic directions in the phase space 
of chaotic systems. Here we review the basic results of ergodic theory, with a 
specific reference to the implications of Oseledets' theorem for the properties of 
the CLVs. We then present a detailed description of a "dynamical" algorithm to 
compute the CLVs and show that it generically converges exponentially in time. 
We also discuss its numerical performance and compare it with other algorithms 
presented in literature. We finally illustrate how CLVs can be used to quantify 
deviations from hyperbolicity with reference to a dissipative system (a chain of 
Henon maps) and a Hamiltonian model (a Fermi-Pasta-Ulam chain). 
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1. Introduction 

Lyapunov exponents (LEs) have been recognized as a central tool for the 
characterization of chaotic dynamics since the late '70s when effective algorithms 
have been independently proposed by Shimada and Nagashima |T] and Benettin et 
al. [2]. In fact, although it had been theoretically clarified that LEs provide the right 
framework for an appropriate (coordinate-independent) characterization of dynamical 
regimes (see the pioneering work by Oseledets it was necessary to wait until an 
effective algorithm was made available to a wide community of scientists, to appreciate 
the usefulness of such a tool in physically relevant models, where it is not obvious 
whether theoretical results apply. It is nowadays widely recognized that LEs help to 
quantify a number of interesting physical properties such as dynamical entropies and 
fractal dimensions [4]. 

Since the LEs can be viewed as the eigenvalues of suitable matrices, one might have 
expected a comparable interest to be devoted to the corresponding Lyapunov vectors. 
Unfortunately, for a long time, the only numerically accessible vectors were those 
ones obtained as a by-product of the Gram-Schmidt orthogonalization procedure (an 
intrinsic step of the algorithms defined in ^il2j), which arc coordinate-dependent. As a 
result, Lyapunov vectors did not attract much interest within the nonlinear-dynamics 
community, even though an objective definition of intrinsic, covariant Lyapunov 
vectorf[3 (CLVs) has been given three decades ago [5 . In fact, until recently, as a result 
of the lack of effective algorithms, only a few papers can be found in the literature 
where the proper CLVs have been determined and used to characterize chaotic states 

The situation has drastically changed a few years ago, when efficient algorithms have 
been developed O [TOj [H] , so that many scientists are now aware that CLVs can offer 
information on the local geometric structure of chaotic attractors, as opposed to LEs, 
which are powerful but global quantities. Several papers appeared, where CLVs have 
been successfully employed to better understand many aspects of chaotic dynamics 
[nillSlliailllllelllllllHIIllEOlEIlEg Some of the relevant questions are touched 
in this Special Issue. 

The goal of this paper is to provide a fairly general introduction to the topic, including 
a description of the various algorithms that have been proposed, but mostly focusing 
on the "dynamical" approach that we believe to be the most effective one |^. In 
particular, we discuss the convergence properties of the CLVs by making use of a 
suitable large-deviation function. An additional point that we explore is the minimal 
angle between the unstable and stable manifold, as it provides direct information 
on the hyperbolicity of the underlying dynamics. In particular, we find that angle 
distribution sastisfies a natural scaling law, both in a dissipative and Hamiltonian 
context. 

More precisely, we start with a discussion of the theoretical context where CLVs 
can be properly defined. This is done in Section [21 which contains a short review 
of the basic results of ergodic theory in the generic context of smooth Riemannian 
manifolds. In Section |3l we then review some implications of Oseledets' theorem 
in the more practical context of JR^ and discrete-time dynamics (to which every 
continuous-time dynamics has to be be reduced in computer simulations) for time- 
reversible systems. In Section|4j we introduce covariant Lyapunov vectors, describe in 

fSometimes they are also refereed as characteristic Lyapunov vectors. 
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detail the algorithm introduced in and discuss its application to non-invertible 
dynamics as well. The convergence properties of the algorithm are discussed in 
Section [Sj with reference to coupled Henon maps and Fermi-Pasta-Ulam chains. In 
Section [S] we briefly review other algorithms that have been proposed in the literature 
and mutually compare them. Section[7]is devoted to a discussion of the angles between 
CLVs, essentially focussing on the angle between the stable and unstable manifold in a 
dissipative (coupled Henon maps) and Hamiltonian (Fermi-Pasta-Ulam chain) model. 
We conclude by briefly mentioning other relevant applications and pending questions. 

2. Basic tools from ergodic theory 

The mathematical foundations of ergodic theory provide a solid basis for the study 
of dynamical systems and, in particular, of their stability properties. For what con- 
cerns the application to various fields like statistical and celestial mechanics, plasma 
and accelerator physics, fluid dynamics, meteorology etc., the main contribution of 
ergodic theory stems from the possibility of translating rigorous mathematical results 
into effective algorithms, that allow for a quantitative analysis. It is worth to point 
out that the tools of ergodic theory apply to both continuous (in time) as well as to 
discrete (in time) dynamical systems, i.e. flows and maps, respectively. The study of 
the Lyapunov characteristic exponents, or simply Lyapunov exponents, is one of the 
main chapters of ergodic theory and goes back to the beginning of the XX-th century, 
when the Russian mathematician Aleksandr M. Lyapunov introduced them as indi- 
cators of the stability of a singular point or of a periodic orbit. As mentioned in the 
introduction, here we aim at providing an overview about covariant Lyapunov vectors 
and their practical use in applications. Accordingly, in this section we summarize the 
main results of ergodic theory that are strictly related to the concept of CLV and to 
the possibility of explicitly computing them 0. 



2.1. Dynamical systems, stability and measures 

Oseledets' theorem is a major achievement in the study of the stability of dynamical 
systems. It proves the existence of LE for generic orbits under quite general conditions 
[3]. In practice, this theorem offers the possibility of extending Lyapunov stability 
analysis from fixed points and periodic orbits, to any trajectory of a dynamical system 
(M, L) , either continuous or discrete in time, defined on a Riemannian manifold M of 
dimension N and equipped with a suitable metric. The application L is the evolution 
operator of the dynamical system, i.e. L*{x) = Xt, where Xt is the image at time t of 
the initial condition x at time t — 0. We assume that at each point x G M we can 
identify a tangent space Tx-M. With reference to a vector ^ e TxM we can define the 
expansion rate 

where the symbol || • || is the norm and DL*{x) : T^Al — > TLt{x)-M indicates the 
evolution in tangent space. The Lyapunov characteristic exponent associated to the 



JFor a more general survey about ergodic theory we suggest the contribution by Lai-Sang Young 
in this special issue 
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vector ^ e T^A^ is defined as 

t-i-oo t ll^l 

provided this limit exists. This quantity essentiahy tells us how an infinitesimal 
perturbation dx of the initial condition x G M along the vector e TxM is 
exponentially expanded (or contracted) at large times 

\\5xt\\ ^ \\Sx\\e^^^'^^\ (3) 

The crucial point here is that for a generic (non periodic) trajectory in A4 at any t the 
tangent application DL^{x) applies to different tangent spaces at different times, and 
the very notion of eigenvalues and eigenvectors of this application is undefined, even 
if the limit ([2|) exists. This notwithstanding, a filtration, i.e. a local decomposition 
into suitable subspaces of T^A^ at any x G still exists as in the case of fixed 
points and periodic orbits. In fact, when the vector ^ is varied in TxM the quantity 
A(^,a:) takes a finite number m < N oi distinct values \\{x) > X2{x) > • • • > X^i^)- 
The filtration of the tangent space into subspaces Si, TxM = Si D S2 ■ • • D Sm, is 
such that by choosing £ Si\ S'i+i, one has A(^,x) — A*(a;). It is important to 
point out that in general one can associate to each A* a multiplicity (or degeneracy) 
gi = dim5i — dimS'i+i. A basis of TxM, obtained by taking in each subspace Si \ Si+i 
a number gi of independent vectors is called a normal basis. It can be easily argued 
that if (/i, /2, • • • , /n) is a generic basis of the tangent space TM and (ei, 62, • • • , ejv) 
is a normal basis one has 

N N 



i=l i=l 



\{e,,x)<y\{U,x) (4) 



where the equality holds only if (/i, /2, • • • , /w) is a normal basis. The ordered 
sequence of the m characteristic exponents A* (a:), each one repeated with its 
multiplicity gt is called the spectrum, Sp{x), of characteristic Lyapunov exponents 
Xi{x) > X2{x) > ••• > AAr(x). In practice, the existence of a filtration amounts 
to establishing the existence of linear subspaces Ei — Si\ Si^i of dimension gi that 
identify the characteristic Lyapunov exponent Xi (x) together with its multiplicity. One 
can easily realize that the extension of such a concept to linear subspaces E C TM 
with dimi? — p < N is straightforward. In fact one can define the Lyapunov 
characteristic exponent of order p 

t^oo t VoP(fi, • • • ,5p) 

provided this limit exists, where VoF(---) is the volume of the p-dimensional 
parallelepiped generated by the tangent vectors in its argument. Let us stress that 
the LE defined in equations ^ and ([S]) are independent of the chosen metric. 
After all of these preliminary considerations one can summarize the very content of 
Oseledets' theorem that applies to a dynamical system {M,ij,,L) that preserves the 
measure /i and can be either continuous or discrete in time, invertible or non-invertible. 
The theorem states that [5] 

For almost all x d M and for any subspace E C TM, such that dimE = 
p < N the limit in Eq. exists and is finite and, in particular, the limit in 
Eq. (0) exists and is finite for any tangent vector ^ £ TM. Moreover the 
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spectrum is a measurable function of x and for any subspace E C Tx-M. it 
exists a normal basis (ei, 62, • • • , e^y) such that 



p 




i=l 



Some interesting consequences of this theorem are worth to be Usted. First of all, 
one can easily realize that LE are constant of the motion generated by the evolution 
operator L, thus implying that for an ergodic system they are almost-everywhere 
constant, Xi{x) = Xi- If the dynamical system is invertible and the measure fi 
is conserved, the sum of the Lyapunov characteristics exponents (or the Lyapunov 
exponent of order N) is zero, X'-'^^%M,x) = J2i=iHei,x) = 0. In particular, 
for a Hamiltonian system and, more generally, for any symplectic diffeomorphism 
(i.e., canonical transformation of the manifold Ai onto itself) the spectrum of LE 
is symmetric, i.e. Xi(x) = — Ajv-i+i(a;) for i = l,...,A^/2, where each individual 
exponent Xpf/2-j{x) — for j = 0, ■ ■ ■ ,nc — 1 ii there are Uc constants of the motion 
(one of which is necessarily the energy of the system). Moreover, for any continuous 
time dynamical system whose support does not reduce to a fixed point, at least one 
of the Lyapunov exponents should vanish. 

Finally, it is worth pointing out one further major consequence of Oseledets' theorem 
that is of primary importance for the aim of this manuscript. We make reference, 
for simplicity, to the case of discrete-time evolution, where the definition of the 
adjoint of the evolution operator in tangent space is straightforwarcJ|: given DL*{x) : 
TxM -> TLt{x)M. its adjoint operator {DL^{x))* : TLt[x)M^ TxM represents the 
time- reversed evolution in tangent space of a reversible dynamical system (A^,/i,L). 
By definition the operator {DL*{x))* DL*{x) is a symmetric, positive-definite linear 
operator on TxM. and, with the same hypotheses of Oseledets' theorem, one can prove 
that the following limit exists almost everywhere 



Moreover, its eigenvectors {di{x)^ • • • , dNix)} are a normal basis in a;, while the cor- 
responding eigenvalues {(5i(x), • • • , 5n{x)} are such that In 5; (a;) — Xi{x). Notice that 
the operator ([T]) describes the evolution in tangent space forward in time for a dura- 
tion t and, then, backward in time for the same lapse of time. Accordingly, the initial 
and the final tangent spaces coincide, so that the difficulties inherent the extension of 
the Lyapunov stability analysis to generic trajectories in M are removed. Notice also 
that, while in ergodic systems the eigenvalues of I?(a;) are almost everywhere constant, 
the corresponding eigenvectors do depend on x. 



2.2. Dynamical systems and entropies 

For the sake of completeness it is worth discussing shortly the crucial contributions by 
Ya. B. Pesin to the mathematical theory of LE. In a paper of 1976 23] he generalized 

§Notice en passant that this is the general situation one has to deal with when numerical estimates 
of Sp{x) and of the corresponding eigenvectors have to be performed. In fact, even continuous— time 
dynamical systems have to be integrated by transforming them into discrete-time algorithms, whose 
reliability is primarily related to the conservation of a suitable measure /i and its symmetries and 
to the choice of appropriate integration time steps that guarantee a sufficient sampling all over the 
manifold JVi 



V{x) = lim [{DL*{x))*DL\x)]-t . 



(7) 



Covariant Lyapunov vectors 



6 



the concept of Lyapunov characteristic exponents for a family of mappings that satisfy 
more general conditions than Lyapunov regularity. In practice, the main result is 
the identification of invariant manifolds associated to nonzero Lyapunov exponents 
in non-uniformly hyperbolic systems: this extends to a larger class of dynamical 
models the study of Lyapunov stability for a measure-preserving dynamics. The 
main achievements by Pesin are contained in |24j . In that manuscript he proved that 
the Kolmogorov-Sinai entropy K{L) of any classical dynamical system (A^,/z,L) can 
be expressed in terms of the LE, through the relation 

K{L)=A f VA,(x)dA*, (8) 
Jm , 

where A is a suitable constant and the sum J^t restricted to the positive Xt . One of 
the important consequences of this result for practical applications is that the entropy 
of a dynamical system (one could say, its degree of unpredictability) has a density 
with respect to the measure fj, and this density is related to the sum of the positive 
LE's. For an ergodic system, since the LE are constant almost-everywhere the above 
expression becomes 

+ 

i^(L)=A^A,; (9) 

i 

In the same paper 24 Pesin discussed how the mechanism of Oseledets' splitting (i.e. 
the existence of a filtration of the tanget-space evolution into invariant submanifolds) 
is related to the hyperbolicity of the dynamical system (AA, fi, L). 



3. The geometrical structure of tangent space 

The existence for invertible dynamics of a coordinate-independent local decomposition 
of the tangent-space evolution into covariant subspaces (again, the mechanism of 
Oseledet's splitting) was discussed by Ruelle in a seminal paper of 1979 [5]. Rather 
than surveying the rigorous mathematical treatment, we have decided to describe 
here, in a language accessible also to non mathematicians, the main consequences of 
Ruelle's contribution. For this purpose we shall adopt notations more familiar to an 
audience of readers interested to applications. In practice, rather than referring to 
an abstract dynamical system, {A4, fi^ L{x)), we shall consider explicitly models of 
physical interest, where M is identified with IR^ and the evolution operator L{x) is 
specified by a set of ordinary differential equations 

xt - f (xO (10) 

where the point x G is now represented by the vector = {x], x^, - ■ ■ , x^} e IR^ 
and the continuous dependence on time is made explicit by the subscript t e IR, so 
that the application f : IR^ IR^ is a continuous-time dynamical rule. In this 
section we limit our analysis to invertible dynamical rules, i.e. we assume that for 
any given position in the phase-space there exists only one backward trajectory. Since 
a continuous-time dynamics can be reduced to a discrete-time one with the help of 
a Poincare section [26], we make reference only to the map formalism, where the 
dynamics reads 

x„ = f(™)(xo) , x„+„=f('")(x„)=f(")(x™)=f(")f(")(xo) (11) 
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where x,„ = {x^, z^j, • • • , x^} € JR^ and the time variable is an integer index m G Z. 
Again, the map f'^™-' is assumed to be invertible, i.e. the inverse map f'^^™^(x,„) exists 
and has only one sohition. 

It is worth pointing out that the assumption of the existence of an invariant measure fi 
is a crucial point that demands much care to be properly formulatecj^. In this discrete 
framework it is convenient to consider the Jacobian matrix at x„ 

J(x„) = ^ , (12) 
whose elements are assumed to be and defined as follows 

[J(x,0].. = (13) 

with Xn^ = [x„]i. The evolution operator in tangent space DL^[x) introduced in 
Eq. ([7]) can be represented by the matrix 

k+n-l 
i—n 

where Mfc^„ : IR^ and the lapse of time t has become the integer k. Note also 

that with a slight abuse of notation we write Mfc,„ instead of Mfe(xn). 
It can be easily shown that Mfc^„ satisfies cocycle properties 

Mfc+,,„ = Mfc,„+,Mi,„ , Mo,„ = I ^ M_fc.„ = (Mfe,„_fc)-i . (15) 

The first relation defines the multiplicative property of Jacobian matrices, while the 
last one defines the inverse, i.e. time-reversed, evolution matrix in the tangent space. 



3.1. Oseledets' matrix and Oseledets' splitting 

It is worth reformulating Oseledets' theorem by taking advantage of the invertible 
dynamics and introducing the concepts of forward and backward Oseledets' matrix, 
denoted by the signs " + " and " — " , respectively. Under the same hypotheses as for 
Oseledets' theorem, one can prove that the following limits exist 

^" = ik In [(Mfe,„)^Mfe,„] (16) 

where T indicates the transpose matri>|j]|. Note that for /c > (fc < 0) the matrix 
(Mfc_„)^Mfc^„ evolves forward (backward) a generic tangent space vector from time 
n to time n + fc, and then backward (forward) in time up to time n. Therefore, the 
forward Oseledets' matrix probes the future dynamics of x„, while the backward 

^In the case of symplectic evolution, ergodicity is usually assumed, in spite of examples where it is 
known that the ergodic hypothesis does not hold (e.g., the Fermi-Pasta-Ulam model). Moreover, even 
uniform hyperbolicity or simply hyperbolicity are not guaranteed on rigorous grounds for many models 
of interest for applications. This notwithstanding, many of the results and tools derived from rigorous 
ergodic theory reveal effective and very useful, when applied to dynamical models that, presumably, 
do not fuUfill all the hypotheses from which they were derived, as captured by the celebrated chaotic 
hypothesis |25]. Anyway, one should always keep in mind that a straightforward application of these 
methods may yield "pathological" outcomes that are just a consequence of singular features emerging 
from the highly complicated nature of fi, provided it does exist. 

I We assume that the scalar product is the one associated with the chosen basis. In order 
to consider a more generic scalar product, one should use the adjoint matrix (M^; „)* instead of 
(Mfe,„)T. 
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Figure 1. (color online) (a) Schematic representation of the action of the 
symmetric matrices M^^Mj; „ on tangent space 7^A4 for finite k. (b) Forward 
(+) and backward (— ) Gram-Schmidt evolution II23I I in tangent space for finite 
k. The oriented black line represent the phase space dynamics mOI I. Blue arrows 
represent the forward dynamics, red ones the time reversed, backward one. 



one probes the past (see Fig. [T]). 

The forward and the backward Oseledets' matrices share the same m < N distinct 
eigenvalues Ai > A2 > ... > A^, each with muhiphcity gi (Yld^iSi — N) which 
coincide with the LEs of the dynamical system pH)) . The orthogonal eigenvectors 

{dn^)± (with i — 1, . . . , N) oi H+ and H^, on the other hand, do differ and are not 
invariant under time reversal. As shown in [5,, one can construct Oseledets' splitting 
fin'' (i = 1, . . . , m) by simply intersecting the so-called Oseledets' subspaces 

f2« = (r«)+n(rW)- (17) 

where 

(rl'')+ = (u«)+®...®(u(™))+ (18) 

and 

(r«)" = (u(^^)"®---®(u«)" (19) 

In both formulae (ul*-*)^ with i = 1,...,to denote the eigenspaces of Oseledets' 
matrices defined in (IT5|) . i.e. they are the (orthogonal) subspaces spanned by the 
eigenvectors of the forward and backward Oseledets' matrices. Note that, as a 
consequence of Oseledets' theorem, the subspaces (rl*'')^ are a filtration of the 
dynamics pUj) or of its time-reversal, that is 

fc-^±oo \k\ \\u\\ 

with the nested subspace structure = (rl^V ^ i^n^)^ D ... D (r^™')+ D 

(ri"+^V EE and = (ri™V D {Tt-'Y D ... D (F^)- 3 (r^"))- ^ 0. 

Oseledets' splitting {fln''}i=i^...^m is generically non-orthogonal, and determines a 
measurable decomposition of the tangent space which is independent of the chosen 
norm, covariant with the dynamics and (obviously) invariant under time reversal, that 
is Mfc.„r2„ = ^fc+ri- The Oseledets' decomposition into such covariant subspaces 
exists for any map (jlip. which is continuous and measurable together with its inverse 
and whose Jacobian matrix exists and is finite in each element. If the Jacobian matrix 
is non-degenerate (m = N), it can be easily shown that also the spectrum of LEs 
is non degenerate and all the subspaces fln^ have dimension one. Conversely, if the 
Jacobian matrix is degenerate, each covariant subspace has a dimension equal to the 
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degeneracy gi of the corresponding LE. 

Finally, as a consequence of Eqs. (|17l20p . it easy to show that a generic tangent space 
vector, u„' e for large values of k grows as 

||M,,„u«||^|lu«|lexp^''=. (21) 
that is, its exponential asymptotic growth rate is the LE A^. 



3.2. Gram-Schmidt vectors and Lyapunov characteristic exponents 

The existence of a well defined limit in Eq. (|16p does not allow one to compute 
neither the spectrum of the Lyapunov characteristic exponents nor the Oseledets' 
subspaces and splitting by directly diagonalizing the symmetric matrices H^, since 
they develop exponentially diverging terms leading to numerical overflows. A way 
out, as noted in Refs. [212], is offered by Eq. The idea is to follow the tangent 
space evolution of p-volumes of increasing dimension p, and to recover the LEs Xi 
from Eq. (|4]). While p- volumes may be defined by an arbitrary basis of p independent 

vectors (gn \gn^', . . . ,gn^), in order to evolve volumes of dimension p > 1, one has to 
take into account the fact that almost every initial vector will converge exponentially 
towards the largest expanding one, 

||Mfe,„gW -Mfc,„g(f)|| ^0 for fc» 1 andVz,j, (22) 

thus quickly making any two vectors i and j numerically undistinguishable. To solve 
this problem it is necessary to periodically orthogonalize the tangent space basis during 
time evolution. If one considers an orthogonal matrix whose vector columns form an 

orthonormal basis in tangent space, G„ = ^gn^|gi' , this can be achieved 

by means of QR decomposition every k timesteps, 

Mfc^nG,, = Gk+,i'R.k,n (23) 

where Gk+n is made of orthonormal vectors at Xfc+„ = f('^^(x„) = f('^+"^(xo) 
and the upper diagonal matrix Rfe^„ contains the information obtained by the 
orthonormalization procedure of the non- orthogonal matrix Gn+k ~ lV[fc.„G„. In 
particular, its diagonal elements [Rfc,n]i,i = 7[.*|, are the local growth rates over a 
time k of the orthogonal vectors G„, that is the growth rates of the axes of inertia of 
p- volumes. In numerical applications, k should not be so large as to produce overflows 
or problems of numerical accuracy in the QR decomposition of G„+fc. In systems 
with an invariant ergodic measure fi on A4, where averages over /i can be replaced by 
time averages, the (ordered) Lyapunov spectrum can be thus computed by taking the 
time-average of the logarithms of these diagonal terms, 

^^-A^^^Ei^^^i:u.- (24) 

The off-diagonal nonzero elements [Rk,n\j,i, on the other hand, are obtained by 
projecting each vector column g^'^j, onto the subspace spanned by {gn]^^} with j < i. 
The orthogonal column vectors of G„j evolving under Eq. (|23p are often referred 
to as forward Gram-Schmidt (GS) vectors (from the well-known orthonormalization 
algorithm). In general, they depend on the initial condition x„ from which the above 
procedure has been initialized, but for m — n sufficiently large they converge to a well 
defined direction which only depends on m. 
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This is the main resuh obtained by Ershov and Potapov |27| : under the general 
hypothesis that the considered dynamical system has an invariant ergodic measure, 
one can show that the orthonormal basis {gn''}i=i....,Ar evolved by Eq. (|23l) converges 
exponentially to the eigenvectors {(dl*'')_}i^i jv of the backward Oseledets' matrix, 

l|g!:ifc-(di^U)-MO for A:-^oo. (25) 
From now on we will always refer to GS vectors as to the "asymptotic" oneJ**l. 

An analogous result holds for the backward Gram-Schmidt vectors of the time reversed 
dynamics, which converge to the same index eigenvectors of the forward Oseledets' 
matrix. Obviously, the (reversely ordered) LEs of the time reversed dynamics 
coincide with the forward ones apart for a sign change 

A. = -A, (26) 

and can be numerically computed by the analogous of Eq. (1^^ . Note that by this 
convention, the fastest growing direction of the time-reversed dynamics is associated 
to the iV-th backward GS vector. 

It is also important to stress that the forward and backward Gram-Schmidt vectors 
are not invariant under time-reversal. Furthermore, they are norm-dependent objects 
- as it follows from the norm-dependent nature of GS orthogonalization - with the 
exception of the first forward and first backward GS vectors (i.e. the ones associated 
to the largest expansion rate in both time directions). 

The fact that the forward GS vectors converge to the eigenvectors of the backward 
Oseledets' matrix could seem odd at first sighl ffjl but it is easy to realize that both 
basis in TnM are only determined by information pertaining to the past dynamics of 
x„, as shown in Fig. [TJ Once again, the same heuristic reasoning applies to backward 
GS vectors and the eigenvectors of the forward Oseledets' matrix. 

4. Covariant Lyapunov vectors 

The possibility to characterize and practically compute a norm-independent and time- 
invariant set of local tangent space vectors associated to the LEs allows one to directly 
probe the expanding and contracting directions of a given dynamical system. This 
development has a number of noteworthy applications beyond the area of ergodic 
theory, ranging from ensemble forecasting to statistical mechanics. It is thus natural 

(i) . 

to consider the covariant Lyapunov vectors v,\ , which are the unitary vectors spanning 
the Oseledet's splitting J7k . 

The knowledge of CLVs allows to identify at each point in phase space a vector 
field that has a natural geometrical interpretation: it is the collection of the tangent 
space varieties associated to the stability properties of the dynamical system ruled by 
the LE's. Accordingly, CLVs are objects of primary importance for the study of a 
dynamical system, since they provide all the information about the local geometrical 
structure of the tangent space. We now discuss them in detail. 

** Ershov and Potapov call the large time limit of the GS vectors a stationary basis, a somewhat 
confusing term since they do vary in time. We therefore prefer to use the term "asymptotic" . 

tfln order to avoid this "inversion" between Oseledets eigenvectors and asymptotic GS vectors, 
some author indeed prefer to name "backward" the GS vectors obtained by moving forward in time, 
and viceversa. However, since we prefer to emphasize the operational way used to compute them, 
i.e. by following the forward dynamics, we will stick with our nomenclature. 
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If dim(rin^) = 1, that is, if is non-degenerate, a single CLV vi*-* associated to 
the LE Xi is uniquely defined (up to a phase factor), while for degenerate LEs 
any nonsingular base formed by gi = dmi{^ln'^) covariant vectors Vn"^ £ $7^1 may 
be arbitrarily considered. In the following, unless otherwise stated, we will unfold 
eventual degeneracies of the Lyapunov Spectrum, and consider N LEs Ai > A2 > 
■ ■ ■ > Xn and their corresponding CLVs {'Vn^}i=i^2,...,N- 

In principle, CLVs can be obtained by: i) a "dynamical" approach, which consists in 
first determining the Lyapunov basis via the GS dynamics (1^^ forward in time and 
and then following the time- reversed evolution in the corresponding subspaces [9]) ; 
ii) a "static" approach which consists in computing the corresponding Oseledets's 
subspaces via p^ . p^ and subsequently intersecting them according to Eq. ((T7| 
[7]. The intersection of large matrices may be susceptible to problems of numerical 
accuracy in high-dimensional spaces. Although some variants of the "static" approach, 
have been recently proposed |101lllj . which improve the performance of such a method, 
we believe that the former "dynamical" approach is more efficient and we thereby 
provide a detailed description in following subsection and discuss its stability in 
Section [5] Anyway, for the sake of completeness, in Section [SJ we briefiy illustrate 
also the static algorithms. 



Dynamical algorithm for computing covariant Lyapunov vectors - Formal aspects 

The dynamical algorithm sketched in ^ makes use of the following simple but powerful 
observation. Suppose we have computed the long-time limit of the forward GS vectors 
(i.e. the asymptotic GS basis) along a certain trajectory. Call these GS vectors gn^ . 
Then, the j-th CLV is contained in the subspace formed by the first j asymptotic GS 
vectors. Moreover, a generic vector u G TR^ evolved backward from time n through 
the time- reversed dynamics within the corresponding Oseledets' subspace (rH^)_ will 
converge asymptotically to the true CLVs Vm'' • The very fact that CLVs are computed 
via a stable dynamical rule and not by intersections of large subspaces ensures the 
numerical stability of this algorithm (see also the next section) . 

Our first observation follows trivially from the definitions (|17m9p . It can be expressed 
as 

v«=Ec(f'^^)gi^'^ (27) 

where the c^f'*"* — (gn''|vi*'') with j < i are the CLVs expansion coefficients at time 
n. Note that CLVs are defined up to an irrelevant sign (i.e. only their orientation 
is defined) and that the first CLV coincides with the first vector of the GS basis by 
construction. 

To demonstrate the second point, we first express the CLVs dynamics in a 
convenient form. Define the matrix whose vector columns are the CLVs, V„ ~ 
^vl^''|vl^^| . . . |vi^-'y Covariant vectors evolve according to 

Mfc,„V„ = V„+fcDfc,„ (28) 

where the diagonal matrix Dfc.„ is composed of the local growth factors 7^, — 

||Mfc_„vi*'' II , that is [Dfc^„]ij- = 5ij7^*|j. For finite fc, the logarithms of these growth 
factors are called finite time Lyapunov exponents (FTLE) 

Af-)=ln7« (29) 
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and their time average obviously coincide with the LEs. 

For the sake of simphcity, we assume that the spectrum is not degenerate, so that all 
individual CLVs may be resolved, and we express Eq. (I27p in matrix form. We write 
V„ = G„C„, where the upper triangular matrix C„ contains the CLVs expansion 
coefficients, [C„]j ,j — Cn'^'^ for j < i. Note that since the CLVs have unit norm by 
definition, the expansion coefficients have to be normalized column by column, 

i 

5^(c(f'^))2 = 1 Vz. (30) 
j=i 

Using Eq. we may rewrite Eq. (1^51) as 

(31) 

which implies 

C„ = Rj,^^C„+fcDfc_„ . (32) 

This is our fundamental equation for backward evolution, showing that the CLVs 
- once expressed in the basis of the GS vectors - can be easily evolved backward by 
multiplying the inverted upper triangular matrices R^^^, where Rfc,n is a "by-product" 
of the forward GS dynamics . 

Now consider a generic nonsingular upper triangular matrix Cn+k which, in the GS 
basis, defines a nonsingular sets of vectors {u^*|j.}i=i,...,jv projected in the Oseledets' 
subspaces, such that Un^j. € (r^*^j,)_. We compare its backward evolution with that 

of true CLVs, u^^*^^ (which in the GS vectors basis we write in the matrix form Cn+k)- 
Evolving both sets of vectors backward according to Eq. ([5^ . we have 

Cri = 'Rf.^Cn+k'Dk^n (33) 

Cn = Rj, ^C„+fcDfc_„ (34) 



which implies 



and thus 



■^n+k 



Dfc,nC„^ — C„+fcDfc_„C„"^ (35) 



C„ ^C„ - D^.^^„C„^j,C„+fcDfc_„ . (36) 
It is easy to verify that 

[C„ C„]ij — — ^[Cj^^j,C„+fe]i,j (37) 

^k,n 

where both sides of Eq. (I37p are upper triangular matrices. Thus, the r.h.s. matrix- 
elements are zero for i > j. On the other hand, for i < j we have 

Iki 

— > for fc — > oo and i < j (38) 

^k.n 

due to the nested structure of the Oseledets' subspaces. For fc ^ 1 one has (up to an 
irrelevant sign) 

[C„^C„]ij 6i^j[C^li^Cn+k]j,j = SijQj (39) 
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where qj are arbitrary factors which solely depend on the angles between the initial 
vectors at time n + k. We finally have, for fc ^ 1, 

[C„],,, «[C„],,,q-i (40) 

that is, under the backward evolution ([^51) almost any nonsingular vector set (once 
it belongs to the Oseledets' subspace) converges backward in time towards the true 
CLV, apart from a trivial scale factor. This proves our second point. 
As a final practical remark, we want to point out that numerical estimates of LE 
and CLV are always performed over long but obviously finite time lapses, whereas 
all theorems hold in the infinite time limit. In the next Section, we show that the 
convergence rate is typically exponential, though the rate may be rather small for 
large system sizes. Accordingly, we do not expect any practical problem. 

4-2. Dynamical algorithm for computing covariant Lyapunov vectors - Numerical 
implementation 

We now discuss the numerical aspects of the dynamical algorithm. The first part 
closely follows the standard algorithm introduced by Benettin et al. to compute LEs 
[2], while in the second part, backward evolution is implemented via Eq. (I28|) . 

(i) Forward transient - Given a generic initial condition xq in phase space and a 
generic set of m < TV orthogonal tangent-space vectors {gp 'ji^i,...,™, we first 
evolve the phase and the tangent space dynamics via Eq. and Eq. 
respectively. This transient should last a number n of timesteps sufficient for the 
phase-space trajectory to converge to the ergodic attractor and for the orthogonal 
vectors to converge to the asymptotic GS vectors. 

(ii) Forward dynamics ~ Once a GS basis is reached, we proceed further evolving 
the reference trajectory x„ and GS vectors by k timesteps at once, recording 
in memory both the local GS vectors Gn+hk and the upper triangular matrices 
R-fc,n+(/i-i)fc, ioT h — 1,2, . . . ,t + to. 

(iii) Backward transient - A generic non-singular upper triangular matrix C„^(j+jj,);; is 
generated and evolved backward via Eq. (1^51) for tQ steps (each moving backward 
by k timesteps). The backward transient length kto should be sufficient to 
converge the tangent space initial conditions close enough to the true CLVs 
expansion coefficients in x„+tfc. Columns of Cn+{t+to)k are kept normalized via 
the condition ((30|). 

(iv) Backward dynamics - Finally, the CLVs expansion coefficients may be further 
evolved backward along the trajectory (always being kept normalized) to sample 
the geometrical structure of the ergodic attractor at points Xn+hk, with h = 
0,1,...,^. CLVs expressed in the phase space coordinates reference can be 
obtained by the GS vectors Gn+hk via Eq. (P7)) . 

Several comments are in order: 

• Sampling frequency - Note that k dictates the sampling frequency of CLVs over 
the phase space, reference trajectory -x.n+hk- It can be adjusted as needed, but 
it should be small enough to avoid numerical overfiows and too large numerical 
errors in the Rfc^m and C„i matrices due to exponential growth and contractions. 
The choice of k thus obviously depends on both the stability properties of the 
chosen dynamical system and on the required sampling frequency. In principle. 
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two different values can be chosen, one (ktr) for the forward and backward 
transients, and a second kgiyn < ktr for the forward and backward dynamics. 

Scalability - The algorithm is fully scalable, and if one just needs to compute the 
first m < N CLVs, that is the m most expanding vectors, no other vectors with 
i > m are needed. Indeed, in order to compute the first m CLVs one only needs 
to consider the N x m reduced GS matrices G„ and the corresponding m x m 
reduced upper triangular matrices Rn. If the dynamics is explicitly invertible, the 
same reasoning could be applied to the computation of the last m CLVs, i.e. the 
m most contracting vectors, by simply writing the dynamical algorithm for the 
time-reversed dynamics (provided that the right trajectory has been generated 
forward in time). 

Degenerate Lyapunov spectra - If the Lyapunov spectrum is degenerate, and some 
covariant subspaces flm have a dimension larger than one, the individual 
vectors spanning such a degenerate subspace have no physical meaning. In this 

(i) 

case, the above algorithm will simply return some arbitrary basis of flm which 
may depend on the initial choice of Cn+(^t+to)k- Fo^ ^ more detailed discussion 
of the structure of the degenerate CLVs associated to zero LE in a quasi-one 
dimensional Hamiltonian systems of hard disks see for instance Ref . [18] . 

Memory issues - All existing algorithms need to store a large amount of 
information. Our dynamical algorithm requires to store the 'Rk,m and Gm 
matrices at periodic sampling points during the forward dynamics, in order to 
be reused during the backward evolution. The matrices Rfe,™ needed to run 
the backward dynamics involve m(m -I- l)/2 floating point numbers at each of 
the t + to sampling points, while the Gm matrices involve mN floating-point 
numbers for t sampling points. In terms of floating-point numbers, the total 
memory requirement is thus 

. + 1" 



Mfnt =tm 



N- 



^^^m(™ + l)_ (41) 



This burden can be sensibly reduced if one is not interested in the spatial structure 
of CLVs but just in their relative angles, since the Gm matrices are not needed 
(the relative angles can be computed directly with reference to the GS vectors 
basis). In this case, the memory requirement is just 

AC ^{t + to) . (42) 

These estimates could nevertheless turn into a large amount of memory for long 
sampling times, large dynamical systems and/or when many CLVs are required, 
quickly exhausting the available fast-access memory. Since accessing disk storage 
every k timesteps could be rather time-consuming (and disk memory in itself 
is a limited resource), we suggest an alternative strategy. Once determined the 
maximum amount of data Mb that can be stored in the fast-access memory, 
divide the total forward integration time k{t + to) into rib blocks of length hb k 
such that Mtot/rib < Mb. After running the usual forward transient in order to 
ensure a proper convergence of the GS vectors at time 0, perform a first forward- 
dynamics run, without storing the 'Rk,m and Gm matrices every k timesteps as 
usual, but just saving the current phase space configuration and the GS vectors 
every hbk timesteps at times nhbk, with n = 0, 1, . . . , nf, — 1, that is, at the 
beginning of each block (this means one does not have to run forward in time 
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the last block). These data can be typically stored on disk, as it will not be 
accessed too frequently. Once the first forward run is concluded, perform a series 
of forward and backward iterations block-by-block, starting from the last block. 
Thus, recover the phase space trajectory and GS vectors at time (ni, — l)hi,k, 
evolve them forward up to time ni,hbk, now storing in the fast access memory 
the R/c,m and Gm matrices every k timesteps, and subsequently evolve backward 
a randomly generated nonsingular initial condition Cnhh^k- Next, pass to the 
(rifc — l)-th block, accessing the data stored at time (rif, — 2)/i(,fc, performing 
the forward evolution up to time {rib — l)hi,k. When coming back, however, 
restart from the coefficient computed in the previous block, C(„j^_i)/ij^fe, so that 
the convergence to the true CLVs is not lost. This procedure is finally repeated 
block by block up to the first one, greatly reducing the amount of memory, at 
the only expense of performing twice the same forward and backward dynamics, 
thus doubling the algorithm computational time. 

• Computational time ~ We give order-of-magnitude estimates for the above 
dynamical algorithm in terms of elementary arithmetic operations (basically 
multiplication/division fioating point operations). Suppose we have a dynamical 
system with N degrees of freedom and we are interested in computing the first 



Forward dynamical evolution requires to run both phase and tangent space 
dynamics for k timesteps and then to perform a single QR decomposition 
(or Gram-Schmidt orthonormalization). Special dynamical systems such as 
the ones with strictly local or globally coupled interactions (in the following: 
easy dynamics) require only 0{N) operations for a single step of phase space 
dynamics and 0{mN) for the tangent space evolution. Unfortunately, generic 
dynamical systems with long-range interactions (in the following: hard dynamics) 
typically require 0{N'^) and 0{mN^) operations for phase-space and tangent- 
space dynamicfJH respectively (the latter with a prefactor < 1). This leads 
to 0{mkN) (for easy dynamics) or 0{m k N"^) (for hard dynamics) operations 
to be performed between consecutive QR decompositions, which is an 0{ni?N) 
algorithm in itself. In particular, the stabilized Gram-Schmidt algorithm requires 
~ 2m'^N operations [28] . 

Backward evolution of the CLVs coefficients via Eq. ((32|) does not require any 
explicit matrix inversion, but can be efficiently computed by back-substitution 
algorithms in ^ operations, while the normalization (1501) requires 0{m?) 

operations. Finally, ~ m^N /2 operations are needed in order to express the CLVs 
in the phase space coordinate basis via Eq. (|27l) . 

Since the convergence of GS vectors and CLVs is exponential (see Section [5] 
below), and in typical applications in systems with many degrees of freedom 
(where computational time may become an issue) one is interested to sample 
the geometrical structure of tangent space over the entire ergodic attractor, one 
usually has t^\ and fcdyn ^ m. For easy dynamics (which nevertheless covers 
the large majority of applications) this implies that the forward evolution is 
dominated by the QR decomposition, while the backward dynamics by the back 
substitution, for a total computational time 



m CLVs. 




(43) 



HWe assume the usage of standard matrix multiplication algorithms. 
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Note that even when the full set of CLVs is computed, m — N , the backward 
evolution is approximately 2.4 times faster than the forward one, or 6 times if one 
is not interested in expressing the CLVs in the phase space coordinate basis. The 
situation is slightly less favorable for hard dynamics, which can be dominated by 
the forward evolution in tangent space, requiring full matrix multiplications every 
time step, for a total time Ttot ^ tkm iV^. 

• Transient length - As mentioned in the previous subsection, the convergence rates 
in both the forward and backward transient rates are exponential and related 
to the difference between consecutive Lyapunov exponents (see also Section [5] 
below). For instance, the convergence to the i-th CLV (z > 1) depends on 
the difference Xi — Xi-i- Dynamical systems with many degrees of freedom are 
typically characterized in the N — )■ oo, by a limit spectrum spectrum \{i/N) that 
is piecewise continuous ^29] (possibly after removing those exponents connected 
to collective and sub extensive modes [HI ISO])- This implies that the difference 
between consecutive exponents scales to zero as 1/iV, so that it is advisable, when 
performing a finite size analysis of such systems, to scale the transient time with 
the number of degrees of freedom, tg — >■ toN . 

4-3. Non-invertible dynamics 

We finally discuss the extension of the concept of CLVs and of the dynamical algorithm 
to non-invertible dynamics. Mathematically, the construction of Oseledets' splitting 
requires the uniquely determined future and past images of any phase space point 
x„ on the attractor. Although it is impossible to meet such a requirement in a non- 
invertible dynamical system (such as, e.g. a chain of logistic maps), the difficulty may 
be circumvented by identifying the past images as those ones that have been actually 
visited during the forward evolution leading to x„. In this way, we can "artificially" 
restore the uniqueness of past images, and be free to apply all the machinery described 
above for invertible dynamics. 

There is of course a price to pay: some points x„ on the ergodic attractor can be 
reached by following different past trajectories. This means that CLVs are not uniquely 
defined, but may depend on the past trajectory followed to reach a given point. There 
is no longer a unique CLV in each point in phase space. This is not really an obstacle 
whenever one, instead of being intersted in reconstructing the local tangent-space 
structure, aims at determining general statistical properties of the CLVs. 

5. Convergence towards the covariant Lyapunov vectors 

We now analyze more closely the convergence of the dynamical algorithm towards the 
CLVs, making use of fiuctuations of finite time Lyapunov exponents. Let us assume 
that the forward evolution has been carried out for a long-enough time to ensure a 
proper convergence to the maximally expanding subspaces (i.e. the GS vectors) along 
a given trajectory f("^(xo) for n — 0,1, ... ,t + to, and that the CLVs vi^' are also 
known at each point x„ for n — 0,1, . . . ,t. If io 3> 1 we can assume these CLVs to be 
perfectly converged. 

Now repeat the backward procedure starting for some randomly chosen vector at time 
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n + k^t + to — n, 

(44) 

and evolve backward up to time n. This procedure, allows to study the convergence 
towards vi*' of the "approximate" CLV ul*' (fc) as a function of the finite backward 
evolution time k on which the vector depends. Let us define Sun^ (fc) = \un^ (fc) — vi* | 
which, being both vectors with unit norm, is just a function of their relative angle 

0^:\k), 

(fc) ^ V2 y^l-cos#(fc) (45) 
(note that 6un\k) ~ 0n^(A:) for small angles). 

Since we are interested in the average convergence rate towards the CLVs, we consider 
the generalized distance 

A(fc,g) = ([(5««(fc)]^)i/' (46) 

where (•) denotes an average over both the ending points x„ and the different initial 
conditions c^^j, in tangent space. Moreover, in order to lighten the notations, we have 
dropped the dependence on the index i (that is kept only when strictly necessary). 
Let us now consider a chain of Henon maps |31) , 

X„+l(£) =a - [xnii)+e{xnii-l)-2Xni£)+Xn{e+l))Y + bXn-li£), (47) 

where i = 1, . . . , iV is a lattice index and periodic boundary conditions are assumed. 
We have studied a set of = 5 maps for a = 1.4, b = 0.3, and e = 0.025. For 
such parameter values there are 5 positive Lyapunov exponents and 5 (quite large) 
negative ones. We focus on the the convergence towards the second CLV following the 
procedure sketched above (the first one is automatically generated during the forward 
iteration) . 

In Fig. [2^ we plot the behaviour of the generalized distance (|46l) as a function of 
the backward convergence time for four different q values. We see that the different 
g-distances all converge to exponentially, thus confirming the quick convergence of 
the method. However, different q-distances converge with different exponential rates. 
In order to describe the process, it is reasonable to conjecture that the evolution 
towards the i-th CLV, as given by Eqs. (|37II38I) . is ultimately controlled by the largest 
ratio between the various subspaces growth factors 7^*^. This is tantamount to 
assuming that all other off-diagonal terms have already converged to zero and at 
some time m < n + k the vector u^^} (fc) lies in the plane defined by the {i — l)-th and 
the i-th CLVs. 

In this plane, the dynamics is controlled by the fiuctuations of the finite time difference 
of the two FTLE, A5 = Xi{k) — Xi_i{k) (we implicitly assume that the spectrum is 
not degenerate in i), which for large k converges to the difference between the two 
corresponding Lyapunov exponents, AA = — Ai_i < 0. Therefore, for large times, 
Ujn tend to rotate towards the i-th CLV and their relative angle 9m.{k) approaches 
zero for a large majority of initial conditions. 
We now consider the restricted generalized distance 

e''(fc,g)~(0™(fc)«) (48) 




Figure 2. (color online) Convergence of the second CLV in a chain of 5 coupled 
Henon maps (periodic boundary conditions) as defined in Eq. II47I I for a = 1.4, 
b = 0.3 and e = 0.025. The uncertainty A(A:, q), defined in Eq. Ij46|l . is plotted in 
panel a for different ij- values. The dashed black line marks exponential decay with 
a rate A2 — Ai. In panel b, the corresponding large deviation function is plotted 
(by following the procedure discussed in the text). Finally, the inset contains the 
corresponding curve q{Xs)- 



where the average has been restricted to small angles only. For small angles we may 
also write 

0„,(fc)~^oe^^^ (49) 
Accordingly, the fluctuations of can be described by a large-deviation function [32] 

P(A5,fc)«e-^(^^)'= (50) 
with S'(AA) = being the minimum of S. Thereby, 

e''{k,q)= dXsPiXs,k)0m{ky = dAie'^I^^"-^'^-')! EEe^(«"= (51) 



where, the upper bound to the integration interval is needed to avoid the unphysical 
divergence of the angle. As long as the maximum of the exponential integrated in 
Eq. (|51l) is reached for a negative A^ (which is itself a function of q), that is, if the 
maximum lies inside the integration interval, L{q)/q coincides with the generalized 
Lyapunov exponent (difference) A^ [33], and 

L{q)^Kq = Xsq~S{Xs) , q^S'CXs). (52) 

It is easy to check that for q 0, Ag/q converges to the Lyapunov exponents difference 
AA (which is the extremal point where S{Xs = AA) = 0). 

On the other hand, if S{Xs) extends to the positive semi-axis, there always exists a 
critical 

qc = S'{0) (53) 

above which, the maximum is attained for A5 = and 

L{q) - -5(0) . (54) 

To resume, we have that L(q) is non positive and monotonically decreasing for 
< g < (?c, from L{0) = to L{qc) = — 'S'(O). For q > q^ finally, we have 
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L{q) = -5(0). 

However, this is not yet the end of the story, since so far we have restricted our 
generalized distance to small angles. In fact, there also exists a family of special initial 
conditions 6*0 that are still close to the {i — l)-th CLV - i.e. the "wrong" direction 
- at time k. They are those for which Oq is extremely close to 7r/2. The fraction of 
such trajectories is exp{Xsk), if < 0, and is of order 1 otherwise. Altogether, their 
contribution Q''(k,q) to the total generalized distance is 

e\k,q)^ dA5e^[^^^-^(^^)l =e^(i)'=. (55) 



We are now in a position to evaluate the two contributions combined, 

A,(fc)= [e'(fc,g) + e'(A:,g)]'/^ . (56) 

If gc < 1 then L{1) — L{qc) and the contribution from Oq — t^/^ never prevails. 
If gc > 1, on the other hand, the large angles contributions prevails for g > 1. Note also 
that if the large deviations are restricted to negative values, qc does not exist (being 
formally infinite) and this latter case applies. To sum up, the generalized distance 
decays exponentially as 

gfeL(i)/g fQj. ^ > min(gc,l) 
gfeA<j/q q ^ min((7c, 1) 

where qc is given by Eq. (|53p . Furthermore, for g — )• 0, the exponential convergence rate 
to thci — th CLV (z > 1) is given by the Lyapunov exponent difference AA^ = A^ — Ai_i, 
while it decreases monotonically as q is increased. 

In Fig.[2)D we plot the large deviation function for Xs = X2 — X1 for the above mentioned 
chain of Hcnon maps. In order to minimize finite-size effects, S{Xs) has been estimated 
by comparing the histograms for two different time lengths ki < k2, 

. logP(A^,fc2)-logP(A^,fci) 

S{Xs) = 7 7 ■ 58 

k2 - ki 

The result for ki = 30 and fc2 = 40 is plotted in Fig. [2j). In the inset we plot the 
derivative of S, which coincides with q. The crossing of q{Xs) with the horizontal axis 
identifies the minimum of S, i.e. the long time average of of A^ that turns out to be 
rather consistent with the direct estimate of AA = A2 — Ai (—0.023 instead of —0.025). 
The crossing with the vertical axis identifies q^ which, in this case is qc = 0.3, so that 
we are in the case gc < 1- 

We are now in the position to compare the direct results plotted in Fig. [5^ with 
the prediction of the multifractal analysis. First of all we note that the exponential 
decay for g = is in perfect agreement with the Lyapunov exponents difference A A 
(see the dashed line). As for the other g-distance, since they are all larger than qc 
their exponential decay rates are expected to be equal to S{0)/q. By computing 
the exponential rate for the largest times, we find that they are consistent with the 
predicted 1/q behavior, although with an estimated prefactor value of S{0) ~ 0.008. 
On the other hand a direct estimate from Fig. [21d yields 5(0) = 0.004: the prefactor 
difference is relatively large, but definitely compatible with the amplitude of the finite- 
size corrections. Notice that since the minimum of S{Xs) occurs for a small A^ value, 
S{0) is quadratically small and this implies that a direct estimate is rather problematic. 
To test our prediction in a different system, wc have also studied a Hamiltonian 
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Figure 3. (color online) Convergence of the 3rd and 9th CLVs in a chain of 5 
FPU oscillators and energy density e = 5. (a) Uncertainty A(fc, q) for different 
g-values and both the 3rd (solid lines) and 9th (dashed line) vectors, (b) Large 
deviation function S for the 3rd CLV. In the inset the corresponding curve ^(A^) 
is shown. Additional technical details are mentioned in the text. 



model, namely a chain of Fermi-Pasta- Ulam (FPU) oscillators [34], again with periodic 
boundary conditions, 

= F[qi+i - Qi) - F{qi - Qi^i) (59) 

where F{x) — x — x^. In this model the energy H — J^'li/'^ + (Qi+i ^ 1"^ + 
{qi+i — 90^/4 is conserved. We have studied a chain with N — 5 oscillators (and 
thus 10 degrees of freedom) and energy density e = H/N — 5. The convergence 
of the 3rd and 9th covariant vectors are reported in Fig. [3^ (see solid and dashed 
lines respectively) . The Q-distance convergence of the two vectors is controlled by the 
differences A3 — A2 and Ag — As, respectively. Given the symmetry of the Lyapunov 
spectrum (A^ = A^v+i-i), the two vectors are expected to converge asymptotically in 
the same way, and this is a another way of checking the correctness of our approach. 
In fact, we see that the three pairs of curves corresponding to different g- values values 
become eventually parallel. The initial difference is due to the fact that the 9th vector 
actually evolves in a higher dimensional space, where additional (and faster) relaxation 
processes are present. 

Furthermore, we have studied the large deviation function for the difference between 
the 3rd and 2nd Lyapunov exponent (by comparing finite-time Lyapunov exponents 
for a time ti = 10 and t2 = 20. The results are plotted in Fig. [SJd. There we see that 
the minimum of S{Xs) occurs for A5 = —0.05 that is quite close to the asymptotic 
value AA = —0.053. Moreover, we have Qc « 0.5, i.e. this case is borderline between 
the two above mentioned classes. Finally, the "asymptotic" exponential rates for the 
curves with q ^ I and g = 2 in Fig.[3K are —0.066 and 0.014, hinting at S{Xs) ~ 0.013, 
a value close to the direct observation 0.011 (from the inset in Fig.jSh), thus confirming 
our theoretical argument. 

To summarize, we have shown that our algorithm converges exponentially to the true 
CLVs. However, when decay rates are computed via finite time ensemble averages, only 
the 0-norm decay rate coincides with the difference between two consecutive Lyapunov 
exponents. Higher norm decay with smaller exponential rates due to fluctuations of 
the FTLEs. The norm-dependence can be explained in the context of large deviation 
theory. 
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6. Static algorithms for covariant Lyapunov vectors 

As anticipated, there exist various algorithms for the computation of CLVs. An 
alternative to the dynamical approach described in Sec. H] is represented by "static" 
algorithms, which do not make use of the intrinsic stability when the backward 
evolution is restricted to suitable subspaced, but rather, determine the CLVs as linear 
combinations of either forward or backward GS vectors. More precisely, at each point 
X along a given trajectory one has (for the sake of simplicity, in this section we omit 
the time index) 

v« = E(g+V^'^)g^'^ = E(g^''V<^^>gi''^ (60) 
where g^'' and g^'' are, respectively, the forward and backword asymptotic GS vectors. 



6.1. Wolfe- Samels on algorithm 

A first static algorithm was introduced by Wolfe and Samelson in Ref. 10 . 
From Eq. (I60p . by using of the identity relation 

N 

E(g^^|g^'^)(g^'VV^) = 5.. (61) 

k=l 

and with the help of simple algebraic manipulations, one obtains (see [TU] for more 
details) 

EE(g+ lg-'^)(g-'^|g+V^''^=0 ^<h (62) 

j=l k=i 

where c^j!''''' = (g^^'' jv^''^) is the CLV coefficient expansion on the forward GS basis. 
Eq. (|6ip can be recast in a matrix form, introducing the h x h square matrix 

k—i 

and the vector 



namely. 



y-1.-(gf M'')) k<h, (64) 



Q('>)yW=o. (65) 

Eqs. ([5^ and hold ioi h > 1 (the first CLV trivially coincides with the first 
forward GS vector) and allow obtaining the expansion coefficients of the h-th CLV as 
the kernel of a matrix computed from the first h forward and h—1 backward asymptotic 
GS vectors (the latter being the stable ones for the time-reversed dynamics). 
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6.2. Kuptsov-Parlitz algorithm 

More recently, Kuptsov and Parlitz have introduced a similar procedure, which 
makes use of LU factorization. 

By using the matrix notation introduced in Section |4l Eq. (|60p can be rewritten as 

V = G+C+=G_C_ (66) 

with the plus and minus indices referring to forward and backward dynamics, 
respectively. Since the forward matrix C+ is upper triangular, while the backward 
one C_ is lower tringular, Eq. (pS)) can be easily recast as an LU factorization [23] 

PC+ = C_ (67) 

where 

P = G^G+ 



Once again, if we are only interested in the j'-th CLV, that is, in the j-th column of 
the matrix only the (/i — 1) x /i upper left corner of matrix P is needed. Since 
C_ is lower triangular and the first j entries of its j-ili column are all zeros, we are 
left with the following system of {j — 1) linear homogeneous equations in j variables 

3 

^[P]fe,, [C+],,, = fc = l,2,...,j-l (69) 

i=l 

which defines the vector j up to a rescaling factor. Not surprisingly, one needs the 
first j forward and (i — 1) backward asymptotic GS vectors to obtain the upper left 
part of P which is needed to compute the (forward) expansion coefhcients of the first 
j CLVs. 



6. 3. Comments 

To conclude, all static methods, be the direct intersection method of Eqs. (I17m9p or 
the two refinements briefly discussed above, require the solutions of certain systems of 
linear homogeneous equations, which in turn depend on the forward and backward GS 
vectors. For direct subspace intersection, the first j forward and the last [N — j + 1) 
backward GS vectors are required to compute the first j CLV, while only the first 
j forward and {j — 1) backward vectors are needed by the refined algorithms - an 
improvement if one is interested in either the first or last CLVs. 

However, in our opinion static approaches suffer from a number of disadvantages. 
First of all, one has to compute both forward and backward GS vectors, and thus 
is forced to perform twice vector orthonormalization at each trajectory point where 
CLVs are needed. Note that, as discussed in Section l4?2l vector orthonormalization 
is the computationally most demanding part of the dynamics for both short ranged 
and globally coupled large dynamical systems. Therefore, it is not a good idea to 
double the number of such orthonormalizations with respect to the plain dynamical 
algorithm of Section 01 

A second concern regards the solution of large systems of linear equations, which 
has to be performed by singular value decomposition (SVD) to attain a satisfactory 
numerical accuracy. However, SVD is more time consuming than back substitution 
by a factor 18, as it requires ~ 6m^ operations for an m x m matrix |28] . 
Finally, for what regards memory requirements, only forward GS vectors need to be 
stored by the static algorithm. While this reduces the memory requirement to about 
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2/3 of what needed by the dynamical algorithm to store both G and R matrices, 
this memory advantage is lost whenever one is only interested in the angles between 
vectors, for which the dynamical algorithm only needs to store the upper triangular 
matrices R and thus finds itself in a better position. 

7. Angles between covariant vectors and subspaces 

Covariant Lyapunov vectors provide direct information on the geometrical structure 
of tangent space. In particular, angles and (near)-tangencies between different CLVs 
or their associated subspaces can be used to characterize the dynamical properties of a 
chaotic dynamical system. In most applications of interest, this amounts to studying 
the distribution of such angles over the (ergodic) attractors. While angles are not 
invariant under a generic coordinate transformation, non-singular transformations 
do preserve zero and non-zero angles. Note also that CLVs corresponding to non- 
degenerate LEs may not become completely parallel along a given trajectory. Since 
they evolve via Eq. ((28| , it is clear that should at any point two vectors be completely 
parallel, they would stay the same along the entire trajectory, contradicting the 
non-degeneracy assumption. However, trajectories can pass arbitrarily close to 
such tangent points, resulting arbitrarily small angles. Therefore, relevant physical 
information has to be encoded in the way the angles probability distributions 
approaches the null angle. 

In Refs. Il6 , angles between CLVs have been studied to show that the tangent space 
of generic spatially extended dissipative systems, such as the Kuramoto-Sivashinsky 
equation or the complex Ginzburg-Landau equation, is split into two decoupled 
subspaces. One comprises a finite number of frequently "entangled" CLVs, or physical 
modes, which carry all the relevant information of the trajectory. A second residual 
set is composed of strongly decaying spurious modes which are transversal to the 
"physical" manifold and themselves organized in mutually transversal subspaces. The 
number of physical modes, which is extensive in the system total degrees of freedom N, 
can be interpreted as the number of effective degrees of freedom needed to faithfully 
describe the chaotic dynamics, leading to the conjecture that the physical modes 
constitute a local approximation of the inertial manifold |35) . 

Another interesting issue regards the degree of non-hyperbolicity of a dynamical 
system. In hyperbolic systems [Ml [37| , there exists a direct sum decomposition of 
the tangent space TxM at each point x into three invariant subspaces 

%:M = ® ® . (70) 

The unstable subspace E" is spanned by the CLVs associated to positive LEs, so that 
each vector u S E!,! is exponentially contracted backward in time. Similarly, the stable 
subspace EJ is spanned by the CLVs associated to negative LEs and any u S EJ is 
exponentially contracted forward in time. Finally, E° is associated to the null LEs 
and their corresponding covariant subspacdsl. In particular, hyperbolicity implies that 
the stable and unstable subspaces are nowhere tangent. Dynamics on a hyperbolic 
attractor is structurally stable, i.e., is insensitive to variations of parameters. It 
manifests strong stochastic properties which allow for detailed theoretical analysis 
and useful results such as the shadowing lemma [38]. Violations of hyperbolicity 
may manisfest themselves either as a change along the attractor of the number of 

§§It is non-empty only for continuous time flows and/or in the presence of conservation laws which 
reduce the attractor dymensionality. 
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stable and unstable directions (unstable dimension variability) or through the presence 
of homoclinic tangencies, i.e. points where the stable and unstable manifolds are 
mutually tangent. 

The fact that the overwhelming majority of dynamical systems of any practical use are 
not hyperbolic may seem a serious drawback, but the "chaotic" hypothesis [25] assumes 
that generic chaotic systems, in spite of violations of hyperbolicity, can be treated as 
essentially hyperbolic ones, allowing one to extend a number of useful results to non- 
hyperbolic chaotic systems. In particular, it implies ergodicity of the attractor and 
the existence of well-defined time averages with a probability distribution satisfying a 
large deviation law. 

While the chaotic hypothesis has an evident empirical success, as testified by the 
positive tests of the fluctuation theorem [311 30] , it is nevertheless of great interest to 
assess the degree of violation of hyperbolicity in a given chaotic system, especially in 
relation to its dynamical properties. In the following, we concentrate on homoclinic 
tangencies. 

7.1. Subspace intersection 

As it has been correctly pointed out in Ref. |13) . in order to compute the angle 
between two linear subspaces it is not sufficient to compute all the angles between 
pairs of vectors taken from two bases spanning the two subspaces. Indeed one has to 
consider angles between arbitrary linear combinations of such vectors. As shown 
in [13], this can be taken into account by SVD. Consider two generic subspaces, 
respectively spanned by mi and m2 different CLVs, mi + < m, with m being 
the tangent space dimension. Suppose TO2 > mi. First of all, we organize the CLVs 
from the two subspaces in two different matrices: Ui (of size to x toi), whose columns 
contain the vectors from Vi and U2 (m x mi), with the vectors from V2. Since angles 
can be equally computed in the phase space coordinate basis as well as in the GS 
basis, we choose the latter and consider the proper columns of the matrix C. 
We have to compute the QR factorization of both matrices, 

Ui=QiRi , U2-Q2R2 (71) 

and compose the toi x m2 matrix Qi. There are mi principal angles ^1'' G [0, 7r/2] 
between the two subspaces, and their cosines are given by the principal values s*-*-* of 
Qf Q2, that is 

sW=cos6lW i = l,2,...,mi (72) 
where we have explicated the time index n. 

In our case we are interested in the intersection between the stable E'* and the unstable 
E" manifolds. Note that since the unstable manifold E" is spanned by the first toi 
CLVs, the matrix Ui is upper triangular, and its corresponding orthogonal matrix Qi 
is just the identity matrix, thus simplifying the calculations of principal angles. In 
particular, we are interested in the minimum angle On — min; 9n^. 

7.2. Numerical examples 

In this section, we present two numerical examples: a Hamiltonian system, the FPU 
chain ([5^ with TV oscillators, and a dissipative one, a chain of N Henon maps (ITT)) . 
In the Hamiltonian case we have toi = TO2 = TV — 2, since there are 4 null LEs in 
one dimensional FPU, corresponding to momentum and energy conservation and the 
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Figure 4. (color online) FPU (at energy density e = 10) and Hcnon chains 
results, (a) FPU Lyapunov spectrum as a function of the rescaled index h = 
(i - 0^5)/iV for different system sizes, N = 16,32,64,128. (b) FPU Minimum 
angle 9 timeseries for N = 64. (c) FPU probability distribution P{0) as a function 
of the rescaled minimum angle for different system sizes, = 16, 32, 64, 128. (d) 
Minimum angle rescaled probability distribution for the Henon chain at sizes 
A'^ = 10 and N = 30. Other parameters as in Fig. |2] 



associated symmetries. In the Henon chain case, one has typicaUy mi — m2 = N, and 
there are no zero LEs. We have investigated the tangent-space geometrical structure 
by measuring the minimum angle between the stable and unstable manifolds along a 
trajectory which samples the system's ergodic measure. 

We start discussing the FPU system at an energy density e = H/N ~ 10, 
which is known to be characterized by a well developed chaotic dynamics. Our 
numerical simulations employed periodic boundary conditions and a McLachlan-Atela 
integration algorithm 141], which is well-suited for Hamiltonian systems, and an 
integration step of At = 0.05 (we have verified the stability of our results versus 
At). We employed a transient (both forward and backward) of to = • 10^ time 
units, and sampled the minimum angle On between the stable and unstable manifold 
along 10^ time units (our sampling and orthonormalization rate is one time unit). 
We considered e = 10 and different system sizes between = 16 and N — 128. 
The rescaled Lyapunov spectra for this system are shown in Fig. |4^, (they have been 
computed both along the forward dynamics with the Benettin et al. algorithm as well 
as from the CLVs expansion rates (|28l) . as a check of CLV numerical convergence). 
The stationary time-series corresponding to the minimum angle On between stable and 
unstable subspaces is plotted in Fig. |4}d, where one can see that it is rather irregular. 
We find therefore convenient to reconstruct the probability distribution P{6) (PDF) of 
the minimal angles, (see Fig. lit). By comparing the distribution for different system 
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sizes, one can conclude that the minimum angle scales as 



as testified by the relatively good data collapse of the various PDFs. This is somehow 
intuitive, considered that the number of principal angles is equal to the stable 
and unstable manifold dimension. Notice that the probability distributions vanish 
algebrically (roughly linearly) as — >■ 0, are characterized by a finite maximum and 
decay exponentially fast for large angles. 

This behavior can be compared with results for a chain of Henon maps. Also for this 
dissipative case, the minimum angle scales as but the shape of the PDF is rather 
different: it starts with a finite value for ^ — )■ and decays esponentially as the angle 
is increased (see Fig. |3Ji). Note also that these latter PDFs are characterized by a 
much larger width. 

One can certainly conclude that the dissipative system is characterized by stronger 
violations of hyperbolicity. It is, nevertheless, important to remark that a vanishing 
PDF when 9^0 does not imply hyperbolicity. In hyperbolic systems, the distribution 
of minimum angles is bounded away from zero 9 . A vanishing P{0) implies that 
trajectori es e xist which pass arbitrarly close to a (zero mesure) set of homoclinic 
tangencie j^^ . The different behavior of P{9) in zero reveals, however, important 
differences in the spatial structure of homoclinic tangencies. In particular, if P{9) ^ 
for ^ << 1, the measure of points characterized by angles (between the stable and 
unstable subspace) smaller than a certain threshold 9o will scale as 6q~^^. The 
implications of this scaling law for the global dynamics are the object of current 
investigation. 



8. Concluding remarks 

Covariant Lyapunov vectors provide an intrinsic decomposition of tangent space 
that is invariant under time reversal and independent of the norm (CLVs coincide 
with the Floquet eigenvectors when computed along periodic orbits |42j). For this 
reason, whenever individual directions need be considered, CLVs have to be chosen 
over other vectors, such as the asymptotic GS vectors or the "singular" vectors [7], 
which represent a finite-time version of the former ones, used in ensemble forecast 
applications and which do depend on the norm. Indeed, it has been shown [5] that 
the orthonormalization procedure (j23p introduces, spurious structures in the individual 
CS vectors which have nothing to do with the underlying dynamics. 
It would be interesting to apply CLVs for a better control of the uncertainty in 
nonlinear models such as those used for weather forecast. While variational data- 
assimilation techniques can benefit from the generic knowledge of the system's unstable 
manifold 03], it is conjectured that the most relevant instabilities which affect large 
scale structures in atmospheric models are not those ones associated to the fastest time 
scales (i.e., the largest LEs), but rather to slower instability modes, characterized by 
longer spatial wavelengths. It logically follows that the most important directions 
for data-assimilation coincide with those associated to small positive LEs, which can 
be individually accessed only by computing CLVs. One may therefore hope that the 

^^Here we assume CLVs in TxM are continuous with respect to x. A sufficient condition for the 
continuity of the asymptotic GS base, which in turns imphes the continuity of CLVs is discussed in 
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ability to access the single unstable directions associated to well-defined timescales 
(the inverse of the LE) is of practical use in optimal forecast methods. 
Covariant Lyapunov vectors have a number of further potential applications in 
dynamical systems theory. They can be used to assess the hyperbolicity of the 
underlying dynamics [9l 113) and, as shown in the previous section, can also contribute 
to shed light on the spatial structure of tangent space. In this respect, it is worth 
mentioning that preliminary simulations performed in the FPU chain at lower energies 
(below the strong stochasticity threshold) reveal stronger violations of hyperbolicity. 
Moreover, CLV have been also employed to characterize the collective dynamics of 
large chaotic systems through their localization properties [15], thus allowing to 
establish a connection between microscopic evolution and the emergence of global 
properties. 

Finally, one puzzling finding in Hamiltonian and symplectic systems regards the 
divergence displayed by the power spectrum of the spatial part of the covariant vectors 
associated to the smallest LEs 0. Further work is required to clarify whether this 
property can be related to the so-called hydrodynamic Lyapunov modes [44 or whether 
this 1// behavior has any relation with actual dynamical properties of these systems. 
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